Commentaires :
Chaque bloc (RNA, protéome, cytokynes) contient un nombre distinct de features.
Les échantillons sont partagés et décrits dans sample_info (avec la variable mode : IR vs IS).
data_t2d <- readRDS("/Users/roxanedesvilles/Documents/Pharmacie P5 DESCARTES/6A - M2 Bioinformatique biologique (UPc) 2025-2026/S1/UE Bioinformatique des omiques avancées/UE 8 - Production et gestion des big data en biologie/tp-multiomics/data_t2d.Rds")
# Vérifier les blocs
lapply(data_t2d, dim)## $RNA
## [1] 15 8933
##
## $protein
## [1] 15 233
##
## $cytokyne
## [1] 15 66
##
## $sample_info
## [1] 15 2
# Nombre d'échantillons et de features
n_samples <- nrow(data_t2d$sample_info)
n_features <- sapply(data_t2d[1:3], ncol)
n_samples; n_features## [1] 15
## RNA protein cytokyne
## 8933 233 66
##
## IR IS
## 7 8
Commentaires :
Bloc RNA : 8933 features.
Bloc protein : 233 features.
Bloc cytokyne : 66 features.
Il y a 15 échantillons au total. 7 sont insulino-résistants (IR), et 8 sont insulino-sensibles (IS).
Commentaire : Pour manipuler plus facilement les données, on transforme les blocs en matrices, puis on utilise tidyverse pour faire un pivot.
# Extraction des matrices
rna_mat <- as.matrix(data_t2d$RNA)
prot_mat <- as.matrix(data_t2d$protein)
cyto_mat <- as.matrix(data_t2d$cytokyne)
# Récupération de l'info classe
classes <- data_t2d$sample_info$mode
samples <- rownames(data_t2d$sample_info)
# Fonction pour convertir en tidy format
tidy_block <- function(mat, block_name){
df <- as.data.frame(mat)
df$Sample <- rownames(df)
df$Class <- classes
df_melt <- melt(df, id.vars=c("Sample","Class"))
df_melt$Block <- block_name
return(df_melt)
}
# Tidy pour chaque bloc
df_rna <- tidy_block(rna_mat, "RNA")
df_prot <- tidy_block(prot_mat, "Protein")
df_cyto <- tidy_block(cyto_mat, "Cytokine")
# Combiner tous les blocs
df_all <- bind_rows(df_rna, df_prot, df_cyto)ggplot(df_all, aes(x=Sample, y=value, fill=Block)) +
geom_boxplot() +
facet_wrap(~Block, scales="free") +
theme_bw() +
theme(axis.text.x=element_text(angle=45, hjust=1)) +
labs(title="Distribution par échantillon", y="Expression")Commentaire : Il faut d’abord vérifier si les données ont été log-transformées et/ou centrées-réduites.
# Calcul CV par gène
cv_rna <- apply(rna_mat, 2, function(x) sd(x, na.rm = TRUE)/mean(x, na.rm = TRUE))
# Résumé
summary(cv_rna)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.03214 0.28737 0.61028 0.74212 1.03213 3.87298 30
# Histogramme
hist(cv_rna, breaks=50, main="CV RNA", xlab="Coefficient de variation", col="green")
abline(v=0.4, col="red", lwd=2)Commentaire : La médiane de la série est d’environ 0.61, avec des valeurs allant de 0.03 à 3.87. On choisit un seuil de 0.4 pour conserver la plupart des gènes.
# Filtrage
rna_keep <- !is.na(cv_rna) & cv_rna >= 0.4
rna_filt <- rna_mat[, rna_keep, drop = FALSE]
# Scaling
rna_scaled <- scale(rna_filt)Les données RNA sont normalisées.
# Écart-type par protéine (feature)
sd_prot <- apply(prot_mat, 2, sd, na.rm = TRUE)
summary(sd_prot)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.09588 0.36629 0.58102 0.77467 1.07512 3.29446
# Histogramme
hist(sd_prot, breaks=50, main="SD Protein", xlab="Standard deviation", col="coral")
abline(v=0.4, col="red", lwd=2)Commentaire : Ici, on choisit de calculer le sd et non le CV. En effet, les données du bloc protein ont déjà été log-transformées (valeurs très négatives). Calculer le CV n’aurait pas de sens ici.
La médiane étant de 0.58 environ (pour des valeurs allant de 0.096 à 3.29), on choisit un seuil de sd de 0.4, qui semble un bon compromis pour filtrer suffisamment de gènes peu variables, tout en restant permissif.
# Filtrage
prot_keep <- !is.na(sd_prot) & sd_prot >= 0.4
prot_filt <- prot_mat[, prot_keep, drop = FALSE]
# Scaling
prot_scaled <- scale(prot_filt)Les données protein sont normalisées.
# Calcul CV par gène
cv_cyto <- apply(cyto_mat, 2, function(x) sd(x, na.rm = TRUE)/mean(x, na.rm = TRUE))
summary(cv_cyto)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02105 0.09890 0.13483 0.22516 0.27289 0.98210
# Histogramme
hist(cv_cyto, breaks=50, main="CV Cytokine", xlab="Coefficient de variation", col="#C8A2C8")
abline(v=0.15, col="red", lwd=2)Commentaire : Les valeurs de CV sont positives et petites (de 0.02 Ã 0.98). On choisit un seuil de 0.15, suffisamment permissif sachant la petite taille de ce bloc dans le dataset.
# Filtrage
cyto_keep <- !is.na(cv_cyto) & cv_cyto >= 0.15
cyto_filt <- cyto_mat[, cyto_keep, drop = FALSE]
# Scaling
cyto_scaled <- scale(cyto_filt)Les données cytokyne sont normalisées.
Visualisation des données après preprocessing
# Supposons que rna_scaled, prot_scaled, cyto_scaled existent après filtrage/scaling
tidy_scaled <- function(mat, block_name){
df <- as.data.frame(mat)
df$Sample <- rownames(df)
df$Class <- classes
df_melt <- melt(df, id.vars=c("Sample","Class"))
df_melt$Block <- block_name
return(df_melt)
}
df_rna_scaled <- tidy_scaled(rna_scaled, "RNA")
df_prot_scaled <- tidy_scaled(prot_scaled, "Protein")
df_cyto_scaled <- tidy_scaled(cyto_scaled, "Cytokine")
df_all_scaled <- bind_rows(df_rna_scaled, df_prot_scaled, df_cyto_scaled)
# Boxplot par bloc après traitement
ggplot(df_all_scaled, aes(x=Block, y=value, fill=Block)) +
geom_boxplot() +
theme_bw() +
labs(title="Distribution par bloc (après filtrage et scaling)", y="Expression (scaled)")## [1] 5844
## [1] 163
## [1] 28
## RNA : 5844 features
## Protein : 163 features
## Cytokine : 28 features
Commentaire : Il reste 5844 features après filtrage du bloc RNA, 163 features après filtrage du bloc protein, et 28 features après filtrage du bloc cytokyne.
# Pour RNA
top_rna <- sort(cv_rna[rna_keep], decreasing=TRUE)
head(top_rna, 10) # 10 gènes les plus variables## ENSG00000124721 ENSG00000105429 ENSG00000130340 ENSG00000118680 ENSG00000134909
## 3.872983 3.872983 3.872983 3.872983 3.872983
## ENSG00000156860 ENSG00000143641 ENSG00000213934 ENSG00000221866 ENSG00000142207
## 3.872983 3.872983 3.872983 3.872983 3.872983
## P05154 Q8WZ74 Q5U077 P11586 V9HW56 Q99963 P13591 P11226
## 3.294457 3.075615 2.654122 2.594447 2.531371 2.473422 2.141749 2.098243
## Q08380 Q4LE79
## 2.063518 2.011576
## ICAM1 GMCSF EGF ENA78 LEPTIN IL17F VEGF RESISTIN
## 0.9821021 0.9620879 0.8744051 0.7949401 0.6861287 0.6427149 0.5400899 0.4145250
## VEGFD MCP1
## 0.3665852 0.3659092
Pour évaluer si les gènes les plus variables au niveau transcriptomique sont traduits, les identifiants ENSEMBL des gènes sélectionnés ont été annotés à l’aide de la base org.Hs.eg.db. Les correspondances SYMBOL et UNIPROT ont été extraites via AnnotationDbi::select().
cv_rna_filtered <- cv_rna[rna_keep]
# 1. Annotation : SYMBOL + UNIPROT pour les gènes RNA
res_annot <- AnnotationDbi::select(
org.Hs.eg.db,
keys = names(cv_rna_filtered),
keytype = "ENSEMBL",
columns = c("SYMBOL", "UNIPROT")
)# 2. Identifier les gènes les plus variables
top_rna_names <- names(sort(cv_rna_filtered, decreasing = TRUE))[1:10]
top_rna_annot <- res_annot %>% filter(ENSEMBL %in% top_rna_names)# 3. Vérifier lesquels sont traduits (ont un UNIPROT)
top_rna_annot$Traduit <- ifelse(is.na(top_rna_annot$UNIPROT), "Non", "Oui")##
## Oui
## 63
# 5. (Optionnel) – Exporter un tableau propre
top_rna_annot <- top_rna_annot %>%
dplyr::arrange(desc(ENSEMBL)) %>%
dplyr::select(ENSEMBL, SYMBOL, UNIPROT, Traduit)
top_rna_annot## ENSEMBL SYMBOL UNIPROT Traduit
## 1 ENSG00000221866 PLXNA4 A0A384NYW6 Oui
## 2 ENSG00000221866 PLXNA4 A4D1N6 Oui
## 3 ENSG00000221866 PLXNA4 E9PAM2 Oui
## 4 ENSG00000221866 PLXNA4 Q6UWC6 Oui
## 5 ENSG00000221866 PLXNA4 Q6ZW89 Oui
## 6 ENSG00000221866 PLXNA4 Q8N969 Oui
## 7 ENSG00000221866 PLXNA4 Q8ND00 Oui
## 8 ENSG00000221866 PLXNA4 Q8NEN3 Oui
## 9 ENSG00000221866 PLXNA4 Q9HCM2 Oui
## 10 ENSG00000221866 PLXNA4 Q9NTD4 Oui
## 11 ENSG00000213934 HBG1 D3GKD8 Oui
## 12 ENSG00000213934 HBG1 D9YZU8 Oui
## 13 ENSG00000213934 HBG1 P02096 Oui
## 14 ENSG00000213934 HBG1 P62027 Oui
## 15 ENSG00000213934 HBG1 P69891 Oui
## 16 ENSG00000213934 HBG1 Q549G1 Oui
## 17 ENSG00000213934 HBG1 Q8TDA1 Oui
## 18 ENSG00000213934 HBG1 Q96FH7 Oui
## 19 ENSG00000156860 FBRS J3KNZ9 Oui
## 20 ENSG00000143641 GALNT2 B7Z6K2 Oui
## 21 ENSG00000143641 GALNT2 A0A1L7NY41 Oui
## 22 ENSG00000143641 GALNT2 A0A1L7NY50 Oui
## 23 ENSG00000143641 GALNT2 A8K1Y3 Oui
## 24 ENSG00000143641 GALNT2 B7Z8V8 Oui
## 25 ENSG00000143641 GALNT2 C5HU00 Oui
## 26 ENSG00000143641 GALNT2 Q10471 Oui
## 27 ENSG00000143641 GALNT2 Q9NPY4 Oui
## 28 ENSG00000142207 URB1 D3DSE5 Oui
## 29 ENSG00000142207 URB1 O60287 Oui
## 30 ENSG00000142207 URB1 Q96NX1 Oui
## 31 ENSG00000142207 URB1 Q9NYQ1 Oui
## 32 ENSG00000134909 ARHGAP32 A0A2X0SFD0 Oui
## 33 ENSG00000134909 ARHGAP32 A7KAX9 Oui
## 34 ENSG00000134909 ARHGAP32 I7H0B0 Oui
## 35 ENSG00000134909 ARHGAP32 O94820 Oui
## 36 ENSG00000134909 ARHGAP32 Q86YL6 Oui
## 37 ENSG00000134909 ARHGAP32 Q8IUG4 Oui
## 38 ENSG00000134909 ARHGAP32 Q9BWG3 Oui
## 39 ENSG00000134909 ARHGAP32 A0A804HK06 Oui
## 40 ENSG00000130340 SNX9 B2RAU5 Oui
## 41 ENSG00000130340 SNX9 Q9BSI7 Oui
## 42 ENSG00000130340 SNX9 Q9BVM1 Oui
## 43 ENSG00000130340 SNX9 Q9UJH6 Oui
## 44 ENSG00000130340 SNX9 Q9UP20 Oui
## 45 ENSG00000130340 SNX9 Q9Y5X1 Oui
## 46 ENSG00000130340 SNX9 A0A7P0T8Z7 Oui
## 47 ENSG00000124721 DNAH8 A0A075B6F3 Oui
## 48 ENSG00000124721 DNAH8 H0Y7V4 Oui
## 49 ENSG00000124721 DNAH8 O00438 Oui
## 50 ENSG00000124721 DNAH8 Q5JYI2 Oui
## 51 ENSG00000124721 DNAH8 Q5T2M3 Oui
## 52 ENSG00000124721 DNAH8 Q5T2M4 Oui
## 53 ENSG00000124721 DNAH8 Q5TG00 Oui
## 54 ENSG00000124721 DNAH8 Q96JB1 Oui
## 55 ENSG00000124721 DNAH8 Q9UEM4 Oui
## 56 ENSG00000118680 MYL12B D3DUH6 Oui
## 57 ENSG00000118680 MYL12B O14950 Oui
## 58 ENSG00000118680 MYL12B Q13182 Oui
## 59 ENSG00000118680 MYL12B Q53HL1 Oui
## 60 ENSG00000118680 MYL12B Q7Z5Z4 Oui
## 61 ENSG00000105429 MEGF8 A8KAY0 Oui
## 62 ENSG00000105429 MEGF8 O75097 Oui
## 63 ENSG00000105429 MEGF8 Q7Z7M0 Oui
traduit_df <- top_rna_annot %>%
count(Traduit)
ggplot(traduit_df, aes(x = Traduit, y = n, fill = Traduit)) +
geom_bar(stat = "identity", width = 0.5) +
geom_text(aes(label = n), vjust = -0.4, size = 5) +
scale_fill_manual(values = c("Oui" = "mediumorchid", "Non" = "grey70")) +
theme_bw() +
labs(
title = "Proportion de gènes RNA traduits (présents dans le protéome)",
x = "",
y = "Nombre de gènes"
) +
theme(
plot.title = element_text(size = 13, face = "bold", hjust = 0.5),
legend.position = "none"
)Commentaires : Après filtrage, tous les gènes RNA les plus variables ont une correspondance dans la base UniProt. ls sont donc potentiellement traduits et présents dans le bloc protéomique.
Cela reflète une bonne cohérence entre les blocs RNA et protéomique, malgré les différences d’amplitude observées lors de la normalisation.
rna <- rna_scaled
pca_rna <- pca(rna, ncomp = 5, center = TRUE, scale = TRUE)
# Individus colorés par classe
plotIndiv(pca_rna, group = data_t2d$sample_info$mode, legend = TRUE, title = "PCA - RNA")# Évaluer le nombre de PC pertinentes
pca_rna$prop_expl_var$X # vecteur numérique de la variance expliquée par chaque PC## PC1 PC2 PC3 PC4 PC5
## 0.19884291 0.12777437 0.09071776 0.08588039 0.07394191
## PC1 PC2 PC3 PC4 PC5
## 0.199 0.128 0.091 0.086 0.074
## PC1 PC2 PC3 PC4 PC5
## 0.1988429 0.3266173 0.4173350 0.5032154 0.5771573
PC1 et PC2 expliquent le mieux la variance des données (32,6% de variance cumulée expliquée).
# Extraire la variance expliquée par chaque PC
var_expl_rna <- pca_rna$prop_expl_var$X # vecteur numérique
# Elbow plot
plot(var_expl_rna, type="b", pch=19, col="steelblue",
xlab="Composante principale", ylab="Variance expliquée",
main="Elbow plot PCA RNA")
abline(h=0.09, col="red", lty=2) # seuil optionnel pour voir les PCs importantes# Récupérer les coordonnées PCA des échantillons
pca_scores_rna <- as.data.frame(pca_rna$variates$X)
pca_scores_rna$sample <- rownames(pca_scores_rna)
pca_scores_rna$class <- data_t2d$sample_info$mode # ta colonne de classes IS/IR
# Plot PCA avec ellipse
library(ggrepel)
pca_scores_rna$outlier <- FALSE
center_IR <- colMeans(pca_scores_rna[pca_scores_rna$class=="IR", c("PC1","PC2")])
dist_IR <- sqrt((pca_scores_rna$PC1 - center_IR[1])^2 + (pca_scores_rna$PC2 - center_IR[2])^2)
pca_scores_rna$outlier[which.max(dist_IR)] <- TRUE
ggplot(pca_scores_rna, aes(x=PC1, y=PC2, color=class)) +
geom_point(size=3) +
geom_point(data=subset(pca_scores_rna, outlier==TRUE),
color="red", size=4, shape=8) + # point rouge étoilé
geom_text_repel(data=subset(pca_scores_rna, outlier==TRUE),
aes(label=sample), color="red") +
stat_ellipse(level=0.95) +
theme_minimal() +
labs(title="PCA RNA avec outlier", x="PC1", y="PC2")Commentaire de l’ellipse :Les deux groupes IS et IR ont une variance similaire sur PC1 et 2. La dispersion des échantillons est comparable. Le groupe IS est plus homogène (l’ellipse de IS est incluse dans celle d’IR). Cette PCA montre que les principales sources de variations ne distinguent pas parfaitement les deux groupes. D’autres blocs pourraient mieux montrer cette séparation.
L’outlier pourrait avoir un profil transcriptionnel un peu différent.
prot <- prot_scaled
pca_prot <- pca(prot, ncomp = 5, center = TRUE, scale = TRUE)
# Individus colorés par classe
plotIndiv(pca_prot, group = data_t2d$sample_info$mode, legend = TRUE, title = "PCA - RNA")# Évaluer le nombre de PC pertinentes
pca_prot$prop_expl_var$X # vecteur numérique de la variance expliquée par chaque PC## PC1 PC2 PC3 PC4 PC5
## 0.22065477 0.15578498 0.12771023 0.08940879 0.08368956
## PC1 PC2 PC3 PC4 PC5
## 0.221 0.156 0.128 0.089 0.084
## PC1 PC2 PC3 PC4 PC5
## 0.2206548 0.3764398 0.5041500 0.5935588 0.6772483
# Extraire la variance expliquée par chaque PC
var_expl_prot <- pca_prot$prop_expl_var$X # vecteur numérique
# Elbow plot
plot(var_expl_prot, type="b", pch=19, col="steelblue",
xlab="Composante principale", ylab="Variance expliquée",
main="Elbow plot PCA RNA")
abline(h=0.155, col="red", lty=2) # seuil optionnel pour voir les PCs importantes# Récupérer les coordonnées PCA des échantillons
pca_scores_prot <- as.data.frame(pca_prot$variates$X)
pca_scores_prot$sample <- rownames(pca_scores_prot)
pca_scores_prot$class <- data_t2d$sample_info$mode # ta colonne de classes IS/IR
# Plot PCA avec ellipse
ggplot(pca_scores_prot, aes(x=PC1, y=PC2, color=class)) +
geom_point(size=3) +
stat_ellipse(level=0.95) + # ellipse 95% pour chaque groupe
theme_minimal() +
labs(title="PCA Protein", x="PC1", y="PC2")Commentaire de l’ellipse :Pour la PCA des protéines, les deux premières composantes principales expliquent une proportion plus élevée de variance que pour les données RNA, ce qui indique que la structure globale des protéines est mieux capturée par les deux premières PC. Les ellipses des groupes IS et IR sont horizontales, mais ont des tailles différentes : le groupe IR présente une dispersion plus large, alors que le groupe IS est plus concentré. Contrairement aux RNA, aucun point n’est un outlier et les chevauchements sont limités, ce qui traduit une séparation plus nette des groupes dans l’espace des deux premières composantes.
cyto <- cyto_scaled
pca_cyto <- pca(cyto, ncomp = 5, center = TRUE, scale = TRUE)
# Individus colorés par classe
plotIndiv(pca_cyto, group = data_t2d$sample_info$mode, legend = TRUE, title = "PCA - RNA")# Évaluer le nombre de PC pertinentes
pca_cyto$prop_expl_var$X # vecteur numérique de la variance expliquée par chaque PC## PC1 PC2 PC3 PC4 PC5
## 0.31421153 0.27689350 0.16856305 0.11865865 0.03686007
## PC1 PC2 PC3 PC4 PC5
## 0.314 0.277 0.169 0.119 0.037
## PC1 PC2 PC3 PC4 PC5
## 0.3142115 0.5911050 0.7596681 0.8783267 0.9151868
# Extraire la variance expliquée par chaque PC
var_expl_cyto <- pca_cyto$prop_expl_var$X # vecteur numérique
# Elbow plot
plot(var_expl_cyto, type="b", pch=19, col="steelblue",
xlab="Composante principale", ylab="Variance expliquée",
main="Elbow plot PCA RNA")# Récupérer les coordonnées PCA des échantillons
pca_scores_cyto <- as.data.frame(pca_cyto$variates$X)
pca_scores_cyto$sample <- rownames(pca_scores_cyto)
pca_scores_cyto$class <- data_t2d$sample_info$mode # ta colonne de classes IS/IR
# Plot PCA avec ellipse
ggplot(pca_scores_cyto, aes(x=PC1, y=PC2, color=class)) +
geom_point(size=3) +
stat_ellipse(level=0.95) + # ellipse 95% pour chaque groupe
theme_minimal() +
labs(title="PCA Cytokine", x="PC1", y="PC2")Commentaire : Pour la PCA des cytokines, les deux premières composantes expliquent moins de variance que pour les protéines, mais restent informatives. Les ellipses des groupes IS et IR se croisent, indiquant un certain chevauchement entre les profils des deux groupes. Les différentes inclinaisons des ellipses IR et IS traduisentt des orientations différentes de la variation dans l’espace des deux premières PC. Aucun point n’apparaît comme outlier, suggérant que toutes les observations sont relativement cohérentes avec la structure globale des groupes.
# Préparer les blocs
# Y = classe IS / IR
Y <- factor(data_t2d$sample_info$mode)
# Les blocs
blocks <- list(
RNA = rna_scaled,
Protein = prot_scaled,
Cytokine = cyto_scaled
)
# sPCA multi-blocs supervisée
spca_res <- block.splsda(
X = blocks,
Y = Y,
ncomp = 2, # on garde 2 composantes
keepX = list(
RNA = c(10, 5), # 10 features sur comp1, 5 sur comp2
Protein = c(10, 5),
Cytokine = c(10, 5)
),
scale = TRUE
)
# Visualiser les échantillons
plotIndiv(spca_res, group = Y, legend = TRUE, ellipse = TRUE, title = "sPCA multi-blocs")# Extraire les variables retenues
selected_vars <- lapply(names(blocks), function(bloc) {
lapply(1:2, function(comp) {
sv <- selectVar(spca_res, block = bloc, comp = comp)
sv[[bloc]]$name # $name pour chaque bloc
})
})
names(selected_vars) <- names(blocks)
# Affichage
selected_vars## $RNA
## $RNA[[1]]
## [1] "ENSG00000145388" "ENSG00000086548" "ENSG00000146858" "ENSG00000136295"
## [5] "ENSG00000100461" "ENSG00000108588" "ENSG00000146556" "ENSG00000165480"
## [9] "ENSG00000185787" "ENSG00000111725"
##
## $RNA[[2]]
## [1] "ENSG00000095564" "ENSG00000125503" "ENSG00000143033" "ENSG00000163607"
## [5] "ENSG00000072832"
##
##
## $Protein
## $Protein[[1]]
## [1] "P01023" "P01019" "Q8IZZ5" "Q92520" "P27918" "P27169" "P02749" "P05154"
## [9] "P36980" "Q96NL6"
##
## $Protein[[2]]
## [1] "P35542" "P02655" "Q9NZP8" "P14151" "Q08380"
##
##
## $Cytokine
## $Cytokine[[1]]
## [1] "VEGFD" "EOTAXIN" "MCP1" "PDGFBB" "EGF" "ICAM1" "ENA78"
## [8] "IL17A" "IL17F" "VEGF"
##
## $Cytokine[[2]]
## [1] "RESISTIN" "RANTES" "LEPTIN" "GMCSF" "IL4"